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This  thesis  investigates  the  feasibility  of  micro¬ 
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minutes.  The  spatial  impulse  responses  for  classical 
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presented  graphically,  along  with  new,  innovative,  spatial 
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Due  to  diffraction  effects,  acoustic  imaging  and  tissue 
characterization  techniques  require  information  about  the 
insonifying  wave  at  the  object's  location.  The  propagation 
of  a  pulsed  ultrasound  wave  with  arbitrary  temporal  and 
spatial  shape  is  not  well  understood,  and  requires 
development  of  computer-aided  analysis  and  simulation 
techniques . 

A  computationally  efficient  method  to  model  ultrasound 
propagation  from  a  planar  source  within  a  medium  has  been 
developed  by  Guyomar  and  Powers  [Refs.  1,2].  The  method  is 
applicable  to  lossless  media,  media  with  a  linear  frequency 
loss  coefficient  such  as  biological  tissue,  and  media  with  a 
quadratic  frequency  loss  coefficient  such  as  liquids  and 
gases.  Relying  on  Fast  Fourier  Transforms  (FFTs)  instead  of 
integral  solutions  that  require  geometric  interpretation 
makes  the  method  very  suitable  for  computer  implementation. 

Models  for  all  three  cases  have  been  previously 
developed,  coded  in  FORTRAN,  and  simulated  on  an  IBM  3033 
mainframe  computer  [Ref.  3].  The  purpose  of  this  thesis  was 
to  investigate  the  feasibility  of  generating  micro-computer 
solutions  of  the  models  within  a  realistic  amount  of  time. 
As  the  average  execution  time  was  thirty  seconds  on  the 
mainframe  for  a  64  x  64  array  size,  a  time  limit  of  thirty 


minutes  was  set  as  a  realistic  goal.  While  the  difference 
between  thirty  seconds  and  thirty  minutes  may  appear 
unacceptable,  rapid  advances  in  micro-computer  technology  is 
narrowing  this  margin  appreciably. 

Since  the  previously  written  FORTRAN  programs  were  un¬ 
available,  the  project  required  starting  from  the  beginning. 
A  first  attempt  was  made  using  a  commercially  available  pro¬ 
gram  called  MATLAB.  HATLAB  is  an  array-oriented  program 
containing  built-in  functions  such  as  Bessel  and  Fast 
Fourier  Transforms,  as  well  as  the  ability  to  program  new 
functions.  Limitations  later  discovered  in  MATLAB  required 
starting  over,  programming  the  entire  project  in  Microsoft 
Fortran  Version  3.31.  The  result  is  a  stand-alone, 
user- interactive  program  allowing  customized  computer  runs 
tailored  to  a  specific  set  of  input  variables.  Only  the 
lossless  and  linear  loss  coefficient  cases  have  been  imple¬ 
mented,  the  third  case  with  quadratic  loss  coefficient  is 
left  for  future  expansion.  The  goal  of  thirty  minutes  was 
obtained  by  code  optimization  and  exploiting  array  symmetry 
where  it  existed. 

The  program  was  tested  and  results  verified  using  clas¬ 
sical  problems  such  as  a  s‘qi»re  or  circular-shaped  piston 
spatial  excitation.  Once  verified,  the  program  was  used  to 
compute  the  spatial  impulse  response  of  new  excitation 
shapes.  Examples  of  these  include  pyramid-shaped  pulsed 
sources  and  pulsed  linear  array  elements. 


The  problem  solved  by  Guyomar  and  Powers  [Refs.  1,2] 
uses  the  geometry  of  Figure  1. 


Given  a  separable  source  of  excitation 

v(x,y,0,f)  =  s(x,y)T{t) ,  (1) 

the  problem  is  to  find  the  acoustic  potential  4>  (x  ,y ,  z  ,t)  at 
an  observation  point  located  a  distance  z  above  the  source 

t 

plane.  A  rigidly  baffled  source  plane  and  linear,  i 

! 
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homogeneous  media  are  assumed.  It  has  been  shown  [Ref.  2] 
that  the  potential  is  given  by 

=  s(x,y)T(t)*x  \  *tg(x,y,z,t) .  (2) 

The  *  indicates  convolution  over  the  variable  appearing 
directly  below  it  and  g(x,y,z,t)  is  the  impulse  response  or 
Green's  function  that  solves  the  wave  equation  and  meets  the 
applicable  boundary  conditions,  as  shown  in  the  block 
diagram  of  Figure  2. 


Figuxa  2.  Block  diagram  of  impulse  response. 


In  the  spatial  domain,  4>(x,y,z,t)  is  also  given  [Refs. 
4-8]  as 

<t>{ x,  y,  z,  t)  =  T(t).;  h(x,  y,  z,  t)  (3) 

where  h(x,y,z,t)  is  the  "spatial  impulse  response"  and  is 
equal  to 

h(x,y,z,t)  =  s( x,y)*x  *vg(x,y,z,t)  (4) 

where  g(x,y,z,t)  is  the  Green's  function  (or  impulse  re¬ 
sponse)  .  In  a  linear  system,  the  relationship  between  the 


spatial  impulse  response,  h(x,y,z,t),  and  the  Green's 
function,  g(x,y,z,t),  is  illustrated  by  the  block  diagram  of 
Figure  3.  Using  the  two-dimensional  spatial  Fourier  trans- 


- 

h(x.y.z.t!  = 
sfc<.yV**glx.y.z.tl 

Figure  3.  Block  diagram  of  spatial  impulse  response. 

form  of  Equation  4,  <£(x,y,z,t)  can  be  written  as 

<Kx>y ,z,t)  =  T(t);  T~l  {sg}  (5) 

where  the  tilde  notation  indicates  the  transform  of  the  spa¬ 
tial  function. 

In  the  spatial  frequency  domain,  the  transform  of  the 
spatial  impulse  response  is 

h(x,y,z,t)  =  T~l  {sg}  .  (6) 

The  convolutions  in  Equations  2  and  4  are  difficult  to 
implement  on  a  computer.  The  technique  presented  by 
Equations  3,  5,  and  6,  and  shown  graphically  in  the  block 

diagram  of  Figure  4,  reduce  two  of  the  convolutions  to 
multiplication. 

> 

i 

i 

I  11 


s(x.y)<5  (t) 


•‘,l4  >  »*!'*•?*%*  \  v  I*  t 


From  Equation  6  we  see  that  if  we  interpret  s  as  the 
angular  spectrum  of  s,  then  Tf  can  be  thought  to  be  the 
"propagation  transfer  function"  associated  with  propagation 


<Mx.y.z.t) 

=  Tit)  *  h(x.y.z.t) 

=  s&.ylTIt)  -#-*•) f  g{x.y.2.t) 
xyt 

Figure  4.  Block  diagram  of  general  solution. 


in  the  medium.  This  propagation  transfer  function  behaves  as 
a  time-varying  spatial  filter  that  increasingly  attenuates 
the  higher  spatial  frequencies  of  the  source  as  time 
increases . 

Three  propagation  models  have  been  developed,  one  for 

♦ 

each  type  of  media.  The  general  technique  followed  for  each 
model  begins  with  the  representative  wave  equation  for  the 
specific  medium  and  finding  the  Green's  function  (or  the  two 
dimensional  spatial  transform  of  the  Green's  function)  which 
solves  the  wave  equation.  Using  this  Green's  function  or 
its  spatial  transform,  the  spatial  impulse  response, 
h(x,y,z,t),  is  found  from  Equation  6.  From  this,  the  acous¬ 
tic  potential  <£(x,y,z,t)  can  be  calculated  for  an  arbitrary 
plane  in  the  positive-z  half  space  using  Equation  5. 


A.  CASE  1:  LOSSLESS  MEDIA 

The  results  of  the  transfer  function  approach  for  loss¬ 
less  media  follow,  as  published  by  Guyomar  and  Powers  [Ref. 
1]. 

The  wave  equation  is  the  Helmholtz  equation 


v2t 


1  d2<f> 
c2  dt* 


=  0. 


The  Green's  function  is 


(7) 


.  .  6{ct  -  R) 


(8) 


where 


R  =  \Jz2  +  y 2  +  z2. 

For  this  problem  the  spatial  itapulse  response  is 

2e(a:,y)  *  -  (K/c)] 


(9) 


y,  0  = 


2  *R 


(10) 


Taking  the  spatial  transform  yields 


9pi  =  [P^rr?)  H{ct  ~  z)  >  (1A) 

where  the  tilde  indicates  the  two-dimensional  spatial 


transform  and 


The  multiplicative  constants  can  be  dropped  since  the  re¬ 
sults  will  be  normalized  to  maximum  values.  The  transfer 
function  for  propagation  in  lossless  media  is  then 


9Pi  =  Jo  (, p'JcH 2  -  z1 )  ff(ct  -  z) ,  (13) 

where  JQ  is  the  Bessel  function  of  the  first  kind  for  order 
zero.  For  programming  considerations  the  term 

z'  =  VcPt2  —  za  (14) 

is  calculated  as  a  unit  and  identified  as  ZPRIME  in  both 
this  paper  and  the  source  code. 

For  a  given  value  of  z,  one  can  calculate  the  value  of 
cfp!  for  each  value  of  time,  then  take  the  inverse  spatial 
transform  of  the  product.  The  result  is  the  spatial  impulse 
response  of  Equation  8. 

B.  CASE  2:  liOSS  COEFFICIENT  LINEAR  IN  FREQUENCY 

For  media  which  exhibits  a  loss  coefficient  that  in¬ 
creases  linearly  with  frequency,  the  results  of  Guyomar  and 
Powers  [Refs.  1,2]  follow. 

The  wave  equation  for  media  with  a  linear  loss  coeffi¬ 
cient  is  the  telegrapher's  equation 


For  this  case  it  is  the  spatial  transform  of  the 
Green's  function  that  is  required.  It  has  been  found  toy  an 


alternate  method  [Ref.  9]  to  be 


'  Jo  [Sp1  -  (A*c*/Q>/c2t2  =1?]  -  z)  for  p  <  Ac/ 2 

Jo  [y/P2  -  (A2cz/4)V'c3«s  -  z2]  e~Ae**/2H(ct  -  z)  for  p  >  Ac/ 2 , 


where  IQ  is  the  Modified  Bessel  function  of  the  first  kind 
of  order  zero.  By  definition,  the  propagation  transfer 
function  is  equal  to  the  transform  of  the  Green's  function, 
so 

<7p2(/x,/v,*,*)  = 

/0  \Vp2  ~  (A2c2/4)  y/c2t2  -  z2]  c~Aj‘:>t/1H{ct  —  z)  for  p  <  Ac/ 2 

Lr _ i  '  (17) 

J0  [\/p2  “  (A2c2/4)  \'c2t2  -  z2j  c~A^tl1H(ct  -  z)  for  p  >  Ac/ 2 . 

The  transform  of  the  spatial  impulse  response  is 

0  = 

*(/*>/yKo  f V P2  -  (A2e2/4)  y/c2t2  —  z2 1  e~Aeitf2H{ct  —  z)  for  p  <  Ac/ 2 
5{fx,fy)Jo  [\/p*  —  (A*c2 /4)  Vc2t2  -  z2]  e~Ac2t/2 H(ct  -  z)  for  p  >  Ac/ 2. 


The  algorithm  for  finding  the  spatial  impulse  response 
is  similar  to  the  lossless  case.  The  transform  is  found  of 
the  driving  function  s(x,y)  and  multiplied  by  gp2  evaluated 
at  each  spatial  frequency.  Taking  the  inverse  transform  of 
this  product  yields  the  required  spatial  impulse  response. 


Ill 


A.  PROGRAMMED  IN  FORTRAN 

Fortran  was  the  language  chosen  for  the  program  for 
three  reasons:  1)  the  availability  of  International 

Mathematical  and  Statistical  Library  (IMSL)  [Ref.  10]  rou¬ 
tines  ,  specifically  those  used  to  calculate  Bessel  functions 
and  one-  or  two-dimensional  Fast  Fourier  Transforms,  2)  the 
portability  of  the  source  code  to  other  brands  of  micro¬ 
computers  and  compilers,  and  3)  familiarity  of  the  source 
code  for  other  students  and  professionals.  Microsoft 
Fortran  Version  3.31  was  the  specific  brand  of  compiler  . 
used. 

B.  ADVANTAGE  OF  SYMMETRY 

Most  of  the  arrays  generated  by  the  program  are  sym¬ 
metrical.  This  was  used  as  an  advantage  in  reducing  the 
number  of  calculations  required  and  therefore  the  execution 
time  of  the  program.  As  operation  of  the  program  is  de¬ 
scribed,  the  areas  where  symmetry  has  been  used  are  identi¬ 
fied.  Figure  5  is  provided  to  aid  in  these  explanations  by 
providing  a  "map"  of  an  array  base.  Figure  5  shows  a  64  x 
64  array  base  situated  on  the  X,Y  plane  and  divided  into 
four  quadrants.  Note  that  the  quadrants  are  unequal  in 
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size,  i.e.,  quadrant  II  includes  rows  1  through  33  and 
columns  1  through  33  for  a  33  x  33  "sub-array"  while  quad- 


Figure  5.  Base  array  configuration. 


rant  I  includes  rows  1  through  33  and  columns  33  through  64 
for  a  32  x  32  "sub-array".  Similarly,  when  the  left  side 
(quadrants  II  and  III)  is  to  be  calculated  as  a  whole,  the 
calculations  cover  all  64  rows  and  columns  1  through  33. 
The  different  sizes  are  required  to  ensure  that  only  column 
33  contains  "peak"  information  and  also  that  element  (33,33) 
is  the  physical  center  of  the  array.  The  importance  of  this 
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and  the  reason  why  it  is  required  will  be  explained  later. 
In-line  code  provides  the  required  "flipping"  actions  to 
expand  one  or  two  quadrants  into  a  full  64  x  64  array  of  in¬ 
formation. 

C.  BLOCK  DIAGRAM 

Figure  6  is  a  block  diagram  illustrating  the  infor¬ 
mation  flow  through  the  program.  The  program  is  explained 
block  by  block  in  the  following  text. 

1.  User  input 

The  program  starts  by  obtaining  information  specific 
for  the  problem  to  be  run.  Using  on-screen  prompts,  the 
user  selects  a  problem  name,  filter  type,  excitation  shape 
and  size,  and  values  for  Z,  RHO,  and  MAXTIME.  If  filter  gp2 
is  selected  for  lossy  media,  the  user  is  further  prompted 
for  a  value  of  the  loss  coefficient  A. 

The  problem  name  assigned  must  be  eight  characters 
or  less  and  conform  to  MS-DOS  filename  specifications.  Once 
entered,  the  program  concatenates  three  different  extensions 
to  the  filename.  This  creates  the  names  of  the  three  output 
files  generated  by  the  program. 

The  first  file  created  is  <f ilename>. IMP,  where 
<filename>  is  the  filename  entered  by  the  user.  The  program 
stores  the  64  x  64  spatial  excitation  array,  s(x,y),  in  the 
default  directory  under  this  name  prior  to  taking  its  two 


di.Mnsi.onal  FFT.  It  is  provided  should  the  user  want  to 
graphically  show  the  spatial  excitation  shape. 

The  second  file  created  is  <filename>.TXT.  This  is 
a  printable  text  file  which  summarizes  the  user's  inputs  for 
a  particular  problem. 


Figure  6.  Program  block  diagram. 


The  third  and  last  file  created  is  <f ilename> . DAT. 
This  file  contains  the  final  64  x  64  array  which  is  the  spa¬ 
tial  impulse  response  for  the  given  problem. 
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All  thres  flits  are  in  ASCII  text  format.  They  can 
be  printed  using  the  DOS  PRINT  command,  put  into  a  word 
processing  program  for  printing,  or  put  into  a  graphics 
program  for  conversion  into  a  picture.  All  figures 
generated  for  this  thesis  were  produced  by  translating  the 
<filename>.IMP  or  <filename>.DAT  file  to  MATLAB  format  using 
the  translation  program  supplied  with  MATLAB.  The  converted 
files  were  then  read  into  MATLAB  and  plotted  using  the  MESH 
command.  As  new,  improved,  graphic  routines  for  micro¬ 
computers  are  introduced,  a  future  improvement  would  be 
incorporation  of  graphic  routines  into  the  program. 

As  the  program  is  written,  fourteen  different  source 
excitation  configurations  are  allowed,  nine  of  which  have 
been  implemented.  The  choice  of  fourteen  different  configu¬ 
rations  was  arbitrary,  providing  a  good  selection  of  choices 
and  ease  of  future  expansion.  Each  of  the  remaining  five 
configurations  can  be  implemented  by  adding  its  name  to  the 
screen  formatting  code  and  writing  the  FORTRAN  code  to  gen¬ 
erate  the  excitation  shape.  The  desired  source  excitation 
shape  is  selected  by  entering  the  selection  number  found 
next  to  the  shape's  description. 

The  variable  RHO  represents  the  radial  distance  in 
the  spatial  frequency  domain  as  shown  in  Equation  12.  It  is 
identified  by  the  Greek  character  p  in  Equations  11  through 
13  and  16  through  18  but  will  be  referred  to  as  RHO  to  match 
the  source  code.  RHO  is  used  to  create  a  1  x  64  vector 
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called  RH01  which  in  turn  is  used  to  generate  a  64  x  64  ar¬ 
ray  called  RH02.  The  vector  RH01  is  unique  since  it  must 
start  approximately  at  -RHO  and  increment  until  it  is  equal 
to  0.0  in  RHOl (33).  From  here  it  continues  to  increase  un¬ 
til  it  equals  +RH0  in  RH01(64).  Since  the  vector  is  of  even 
length,  it  cannot  be  formed  simply  by  dividing  the  entire 
range  of  RHO  (-RH0  to  +RHO)  by  64  and  using  this  value  as 
the  increment.  To  do  so  would  create  equal  values  in 
RHOl (32)  and  RH01(33)  which  defeats  the  requirement  that 
RHOl (33)  contain  singular  information,  in  this  case  the 
value  0.0. 

The  largest  value  of  time,  in  seconds,  for  the  time 
dimension  is  represented  by  the  variable  MAXTIME.  For  each 
problem,  time  starts  at  the  instant  that  the  wave  front 
first  arrives  at  the  observation  plane.  This  is  when  time  - 
Z/C,  where  C  is  the  sound  velocity  (1500  meters/second) ,  and 
increments  linearly  until  it  equals  MAXTIME.  To  prevent  a 
run-time  error,  the  user  is  prompted  for  a  value  of  MAXTIME 
which  is  greater  than  Z/C  seconds. 

The  height  of  the  observation  plane  above  the  X,Y 
plane  is  represented  by  Z  (in  meters)  .  It  is  typically  a 
small  number,  several  hundredths  of  a  meter,  or  can  be  set 
to  zero.  A  value  of  0.1  meters  (10  cm)  was  used  for 
debugging  the  program  and  producing  the  results  found  in  the 
Numerical  Simulations  section. 


At  the  completion  of  user  input,  the  program  starts 
to  generate  the  required  vectors  and  arrays.  The  first  of 
these,  the  TIME  vector,  is  initialized  and  then  used  to 
generate  the  ZPRIME  vector.  Both  vectors  are  1  x  64  in  size 
and  represent  a  series  of  time  steps  starting  at  time  =  Z/C 
and  increasing  to  MAXTIME. 

The  program  contains  a  variable  called  STEP.  Its 

default  value  of  4  (which  can  be  changed)  sets  rows  1 

through  4  of  the  final  RESULT  array  to  zero  after  all 

calculations  have  finished.  This  simulates  the  step 

function  H(ct-z) .  Since  any  information  contained  in  the 

first  four  rows,  is  lost,  these  values  need  not  be 

calculated.  For  this  reason,  TIME  starts  at  TIME(STEP+l) , 

or  TIME (5) ,  with  a  value  of  Z/C  seconds,  and  increments 

linearly  until  it  reaches  MAXTIME  in  TIME (64).  The  smooth 

linear  progression  of  time  from  Z/C  to  MAXTIME  is  obtained 

using  the  following  code  segment: 

INC  ■  (MAXTIME  -  TS) /REAL (NROWS -STEP- 1) 

DO  XX  I  *  (STEP+1) , NROWS 

TIME (I)  =  TS  +  (I-STEP-1)  *  INC 
CONTINUE 

where  TS  =  Time  Start  =  Z/C,  NROWS  =  64,  and  STEP  =  4. 

As  an  example  let  MAXTIME  =  1.0E-4  seconds  and  Z  = 
0.1  meters  as  determined  by  the  user.  The  quantity  TS 
equals  6.666667E-5  seconds  and  this  value  is  stored  in 
TIME (5)  .  Each  element  of  TIME  is  increased  by  INC  = 
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5.649718E-7  until  it  equals  MAXTIME  in  TIME (64) .  The 
elapsed  time  for  this  particular  simulation  is  MAXTIME  -  TS 


-  3 . 333333E-5  seconds. 

The  ZPRIME  vector  is  created  using  values  found  in 
the  TIME  vector  as  shown  in  Equation  14 .  Each  value  of  TIME 
is  squared  and  multiplied  by  C2.  Subtracting  Z2  from  this 
product  and  taking  the  square  root  produces  a  ZPRIME  ele¬ 
ment.  Continuing  with  the  previous  example,  two  sample 
ZPRIME  elements  are 

ZPRIME (6)  -  SQRT(TIME(6) 2  *  C2  -  Z2)  -  1.304644E-2 
ZPRIME (7)  -  SQRT(TIME (7) 2  *  C2  -  Z2)  -  1.848934E-2 
The  algorithm  loops  until  ZPRIME (STEP+i)  through  ZPRIME (64) 
have  been  calculated. 

3.  generation  of  BH02  array 

The  next  significant  step  is  formation  of  the  RH02 
array.  First  the  RHOl  vector  is  created  by  calculating  an 
increment  and  base  value  using  the  code 
INC  -  REAL (RHO) /REAL (NCOLS/2-1) 

BASE  -  REAL ( -RHOl)  -  INC 

Assuming  a  value  of  200  for  RHO,  INC  =  6.451613  and  BASE  = 
-206.451613.  This  results  in  a  RHOl  vector  starting  with 
-206.451613  in  RH01(1),  with  each  successive  element  incre¬ 
mented  by  INC.  This  method  ensures  a  value  of  0.0  in 
RHOl (33) ,  -200.0  (-RH0)  in  RH01(2),  and  +200.0  (+RHO)  in 

RHOl (64 ) .  As  with  the  TIME  vector,  simply  setting  the  RHOl 
vector  to  range  from  -RHO  in  RHOl(l)  to  +RHO  in  RHOl  (64) 


would  result  in  identical  values  in  RH01(32)  and  RH01(33). 
This  would  not  meet  the  requirement  that  RH01(33)  contain 
singular  information. 


After  completing  the  RH01  vector,  its  values  are 
used  to  generate  the  RH02  array.  First  the  RHOl  vector 
values  are  placed  in  row  33  of  RH02  and  its  transpose  in 
column  33.  The  program  then  fills  in  the  blank  array 
elements  in  quadrants  II  and  III  using  the  Pythagorean 
theorem  with  the  values  found  in  the  respective  row  33  and 
column  33  positions.  Continuing  with  the  example  using  RHO 
*  200,  the  value  -206.451613  would  be  found  in  elements 
RHO  2 (33,1)  and  RH02(1,33)  (column  major  format).  Using  the 
Pythagorean  theorem,  element  RH02(1,1)  would  receive  the 
value  291.966671.  This  represents  the  value  of  RHO  as 
calculated  radially  from  the  central  positirn  of 
RH02 (33,33) . 

Since  the  RH02  array  is  symmetrical,  only  quadrants 
II  and  III  (left  side)  must  be  calculated.  Although  the 
time  saved  here  is  small,  the  algorithm  using  these  values 
later  in  the  program  needs  only  the  values  found  in  the  left 
half.  Therefore  the  right  side  values  are  not  calculated. 

4.  Spatial  excitation  generation 

Generating  the  spatial  excitation  array  is  the  next 
step.  The  first  five  selections  represent  single  piston 
spatial  excitations  centered  about  position  (33,33)  on  a  64 
x  64  base  array  called  PULSE.  The  five  shapes  available  are 
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the  square  piston ,  circle  piston,  gaussian,  pyramid,  and 
truncated-pyramid.  All  five  shapes  require  the  user  to 


specify  the  width  (or  diameter)  of  the  excitation.  The 
value  of  width  must  be  an  odd  integer  to  ensure  symmetry 
about  PULSE (33 , 33) ,  and  less  than  or  equal  to  63  since  the 
base  array  size  is  64  x  64.  Small  values  between  9  and  15 
provide  excellent  results. 

The  square  and  circle  selections  create  a  simple 
spatial  excitation  with  amplitude  *  1.0.  A  resolution 

problem  exists  with  the  circular  excitation.  Simply,  it  is 
impossible  to  represent  a  circular  shape  using  a  small 
number  of  square  elements.  The  larger  the  diameter,  the 
smoother  the  outer  surface  becomes,  but  with  a  64  x  64  array 
size  the  range  of  available  diameters  is  limited.  The  only 
solution  is  an  increase  in  array  size  to  128  x  128  r 
larger.  This  is  an  excellent  area  for  future  expansion  as 
micro-computers  become  faster  and  less  memory  limited. 

The  Gaussian  selection  creates  a  circular-shaped 
excitation  with  amplitude  following  a  gaussian  distribution. 
The  spatial  excitation  is  truncated  at  the  user  specified 
diameter,  and  its  width  measured  at  the  1/e  point.  This 
excitation,  as  well  as  the  pyramid  and  truncated-pyramid, 
represent  non-piston  sources  with  spatial  variations  in 
amplitude. 

The  pyramid  creates  a  four-sided  tapered  spatial 
excitation  with  a  peak  amplitude  of  1.0  at  the  center  and 


tapering  off  until  it  equals  0.0  at  the  base.  The  size  of 
the  base  is  determined  by  the  input  value  of  width.  For 
example,  a  pyramid  excitation  of  width  ll  would  have  its 
center  at  PULSE (33, 33)  and  occupy  the  region  bordered  by 
rows  29  through  37  and  columns  29  through  37.  The  amplitude 
would  decrease  linearly  and  equal  zero  at  the  outer  edges  of 
the  above  region.  The  truncated-pyramid  is  similar  except 
that  the  amplitude  tapers  off  at  a  slower  rate  until  it 
equals  0.0  at  the  array's  outer  edges.  At  the  specified 
value  of  width,  the  pyramid  is  truncated,  resulting  in  a 
"square  house"  shape  with  tapered  roof  and  square  sides. 

The  remaining  selections  all  represent  linear  arrays 
of  five  elements  each.  An  individual  element  size  of  5  x  5 
located  *  on  10  unit  wide  center  points  was  selected  for 
convenience.  Shape*  selections  for  the  individual 
excitations,  are  square,  circular,  and  gaussian.  Two 
additional  arrays  also  use  5  square  elements,  but  the 
amplitude  of  the  elements  are  stepped  or  parabolic  in  shape. 
The  stepped  array  has  a  center  element  of  amplitude  =  1.0, 
two  adjacent  elements  of  amplitude  =  2/3,  and  two  outer 
elements  of  amplitude  =  1/3.  The  five  elements  thus  create 
a  stepped  appearance:  1/3,  2/3,  1.0,  2/3,  1/3.  The 
parabolic  array  is  similar  to  the  stepped  array,  except  the 
amplitudes  follow  a  parabolic  shape. 

After  the  spatial  excitation  is  generated,  the  two 
dimensional  Fast  Fourier  Transform  is  taken  using  the  IMSL 
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routines  FPT3D  and  FFTCC.  The  routine  FFT3D  computes  the 
FFT  of  a  complex-valued  two-dimensional  array  by  passing 
first  the  columns  and  then  the  rows  as  individual  vectors  to 
the  routine  FFTCC  which  computes  the  FFT  of  the  vector  argu¬ 
ment.  The  array  returned  by  FFT3D  is  complex  and  has  its 
peak  information  located  at  the  outer  four  corners  of  the 
array,  with  the  largest  value  in  PULSE (1,1).  The  array  is 
shifted  to  bring  this  peak  information  to  the  center  using 
the  subroutine  SHIFT.  This  subroutine  swaps  quadrant  I  with 
quadrant  III  and  quadrant  II  with  quadrant  IV,  thus  trans¬ 
ferring  the  largest  value  from  element  PULSE (1,1)  to  element 
PULSE (33, 33) .  This  is  why  it  is  so  important  to  arrange  all 
arrays  so  their  central  value  is  in  element  (33,33). 

5 .  Filter  qpi  lor  loaglggg  media 

From  this  point  the  program  branches  based  on  the 
filter  type  selected.  The  algorithm  for  the  lossless  case, 
using  filter  gp^,  is  to  take  a  value  of  ZPRIME,  multiply  the 
RH02  array  by  this  value  and  pass  each  product  as  the 
argument  to  the  IMSL  subroutine  MMBSJN,  where  MMBSJN 
calculates  the  Bessel  function  of  the  first  kind  of  order 
zero  for  a  real  number  argument.  The  result  returned  from 
the  Bessel  function  is  multiplied  by  the  respective  element 
in  array  PULSE  and  stored  in  array  WORK.  The  next  value  of 
ZPRIME  is  selected  and  the  calculations  again  performed. 


As  this  portion  of  the  program  is  the  most  time 
consuming,  two  specific  actions  were  taken  to  optimize  the 


cod*.  As  previously  explained,  the  number  of  rows 
identified  by  STEP  are  zeroed  in  the  final  array.  Because 
of  this,  calculations  start  at  row  (STEP+l) ,  thus  saving  the 
calculation  of  4  rows  for  each  iteration.  Also,  since  the 
arrays  KH02  and  PULSE  are  symmetrical  in  all  four  quadrants 
about  element  (33,33),  only  quadrant  II  needs  to  be 
calculated.  With  the  default  64  x  64  array  size  and  STEP  * 
4,  only  59  (64-STEP-l)  subarrays  of  size  33  x  33  must  be 
calculated  for  a  total  of  64,251  array  elements.  This 
compares  favorably,  with  the  262,144  array  elements  that 
would  require  calculation  if  the  entire  64  x  64  array  was 
calculated  64  times.  It  is  important  to  note  that  the 
single  most  time-consuming  process  is  the  calculation  of  the 
Bessel  function. 

After  each  33  x  33  subarray  is  calculated,  inline 
code  is  used  to  fill  the  entire  WORK  array.  Quadrant  II  is 
flipped  to  fill  quadrant  III,  then  the  entire  left  side  is 
flipped  to  fill  the  right  side. 

At  this  point  the  inverse  two-dimensional  Fast 
Fourier  Transform  of  the  WORK  array  needs  to  be  taken.  But 
instead  of  passing  the  entire  array  as  an  argument  to  a  two 
dimensional  inverse  FFT  subroutine,  the  array  is  processed 
one  column  at  a  time.  Using  the  IMSL  subroutine  FFT2C, 
which  computes  the  inverse  FFT  (or  FFT)  of  a  complex-valued 
vector,  each  of  the  64  columns  are  passed  individually  as  a 
vector.  After  all  the  columns  are  transformed,  only  one 


row,  row  33,  is  passed  as  an  argument.  Assuming  a  symmetri¬ 
cal  transducer,  row  33  contains  information  from  the  central 
line  across  one  axis  of  the  observation  plane.  The  inverse 
transform  of  this  row  represents  the  field  values  along  this 
line  for  the  specific  value  of  ZPRIME,  and  it  is  this 
information  that  is  used  to  sequentially  build  the  final 
array  called  RESULT.  The  59  iterations  completed  by  the 
above  algorithm  will  each  supply  one  row,  corresponding  to 
an  instant  of  time,  in  array  RESULT.  Specifically,  with  the' 
default  value  of  STEP,  rows  5  (STEP+1)  through  64  (NROWS) 
are  filled. 

After  completing  row  64 ,  the  number  of  rows 
identified  by  STEP  are  set  to  zero  as  they  do  not  contain 
any  information.  The  entire  RESULT  array  is  then  saved  to 
disk,  and  the  program  terminates  normally. 

6 .  Filter  ap2  for  lossy  media 

In  the  second  case,  using  filter  gp2  for  media  with 
linear  loss  coefficient,  the  algorithm  has  two  additional 
steps.  The  program  implements  Equation  17  and  therefore 
must  test  each  argument  for  RHO  greater  than  AC/2.  For 
RHO  >  AC/2,  the  argument  is  passed  to  the  IMSL  Bessel 
function  subroutine  MMBSJN  as  in  the  lossless  case.  But,  if 
RHO  25  AC/2,  the  argument  is  instead  passed  to  the  IMSL 
subroutine  MMBSIN  which  calculates  the  modified  Bessel 
function  of  the  first  kind  of  order  zero.  The  result 


returned  from  either  Bessel  function  is  multiplied  by  the 


respective  element  of  PULSE.  Then,  as  shown  in  Equation  17, 
the  product  is  multiplied  by  an  exponential  term  prior  to 
storing  in  array  WORK.  This  is  the  second  difference  in  the 
lossy  case. 

The  algorithm  loops,  as  in  the  lossless  case,  to 
fill  the  entire  RESULT  array.  After  saving  the  RESULT  array 
on  disk,  the  program  terminates  normally. 
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After  careful  debugging  and  testing  on  a  point  by  point 
basis,  the  program  was  used  to  calculate  the  spatial  Impulse 
response  for  various  source  excitations.  For  each  excita¬ 
tion  shape,  both  the  lossless  and  lossy  case  were  calculated 
and  the  results  plotted.  The  plots  show  the  amplitude  of 
the  wave  plotted  against  cross-direction  and  time  at  an  ob¬ 
servation  plane  located  above  the  X,Y  plane.  In  all  calcu¬ 
lations  the  height  of  the  observation  plane,  Z,  is  10  cm 


above  the  source  plane.  The  value  of  A  used  for  the  lossy 
cases  was  8.0E-3.  This  value  was  selected  as  a  compromise 
and  arrived  at  using  trial  and  error.  Smaller  values  of  A 
caused  only  minor  variations  in '  shape  while  larger  values 
caused  gross  distortion.  All  problems,  unless  noted  other¬ 
wise,  are  calculated  for  a  MAXTIME  of  1.5E-4  seconds.  Thus, 
with  a  Z  value  of  10  cm,  each  result  runs  from  a  starting 
time  »  Z/C  *  6.666667E-5  seconds  to  1.5E-4  seconds,  a  dura¬ 
tion  of  8 . 333333E-5  seconds. 

As  previously  explained,  all  plots  were  obtained  using 
the  MATLAB  MESH  command.  Two  limitations  exist  with  the 
MESH  routine  that  require  further  explanation.  The  first 
limitation  is  the  way  MATLAB  sizes  the  plot.  MATLAB  uses 
input  array  values  to  fill  a  window  of  pre-determined  size. 
This  fixes  the  height  of  the  plot  so  an  excitation  with 
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amplitude  equal  to  1.0  will  plot  identically  to  an 
excitation  with  amplitude  other  than  1.0.  The  second 
limitation  builds  on  the  first:  MATLAB  does  not  show 
numerically  scaled  axes  when  plotting  an  array  of  data.  If 
it  did,  the  different  amplitudes  for  identically  shaped 
plots  would  be  apparent. 

An  alternative  way  to  obtain  numerical  information  was 
found.  Using  the  PLOT  command,  MATLAB  has  the  ability  to 
plot  a  vector  in  the  X,Y  domain  with  numerically  scaled 
axes.  Putting  an  array  of  data  into  PLOT  creates  a  side 
view,  plotting  successive  "slices"  of  the  array  with  time 
increasing  in  the  positive  X  direction.  Using  this  tech¬ 
nique,  amplitude  values  could  be  obtained  for  a  given  plot. 
Since  shape  is  the  primary  interest,  amplitude  information 
obtained  in  this  manner  'was  adequate.  Three  plotted  arrays 
(Figures  9,  11,  and  15)  are  included  for  illustration. 

Execution  times  were  all  less  than  the  stated  goal  of 
thirty  minutes,  with  an  average  of  twenty-five  minutes. 
These  times  were  obtained  using  an  AT&T  6300  micro-computer 
operating  at  8  MHz.  An  8087  numerical  co-processor  chip  was 
installed  and  used  to  increase  performance.  Similar 
performance  can  be  obtained  with  fast  IBM  XT  compatibles  or 
AT  class  micro-computers.  The  next  generation  of 
micro-computers  utilizing  the  Intel  80386  microprocessor  are 
the  logical  choice  for  continued  program  development. 


Figure  7.  Square  piston  spatial  excitation. 


Figure  7  shows  graphically  the  source  excitation  shape 
for  a  square  transducer  of  width  11  and  amplitude  1.0.  The 
excitation  is  centered  over  the  base  array's  central  ele¬ 
ment,  PULSE(33 , 33)  .  In  Figure  8  the  calculated  spatial  im¬ 
pulse  response  for  the  square  transducer  in  lossless  media 
is  shown,  as  is  the  array's  orientation  with  time  and  space. 
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This  orientation  of  time  and  space  applies  to  all  spatial 
impulse  response  plots.  The  origin  of  the  time  axis  starts 
at  Z/C  where  the  potential  is  known  to  be  a  replica  of  the 


excitation  source  shape  for  lossless  media.  As  time 
progresses,  the  potential  is  reduced  until  it  forms  two  out¬ 
ward  traveling  "tails".  As  expected,  the  "tails"  are  square 


Figure  8.  Field  for  square  transducer  (lossless  case) . 

to  match  the  spatial  excitation  shape,  and  slowly  attenuate 
with  time  in  the  lossless  case.  Numerical  information  can 
be  obtained  from  Figure  9  which  is  the  side  view  of  the 
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plot  in  Figure  8.  The  stop  function  is  clearly  visible  as 
is  the  jump  in  potential  to  1.0  in  row  5.  From  this  point 
the  potential  peaks  at  a  value  slightly  greater  than  1.0  and 
then  drops  rapidly.  At  infinity  the  tails  should  reduce  to 
zero. 


REL/  j  TIME 

Figure  9.  Side  view  for  lossless  case. 


In  Figure  10  the  same  excitation  has  been  analyzed  for 
the  lossy  case.  Comparing  the  shapes  in  Figures  8  and  10 
reveals  two  distinct  differences.  The  "tails"  are  no  longer 
square  but  have  taken  on  a  triangular  appearance.  Larger 
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values  of  A  will  produce  a  sharper  triangular  shape.  The 
second  major  difference  concerns  the  region  between  the 
tails.  In  the  lossless  case  this  region  was  at  zero  poten¬ 
tial.  In  the  lossy  case,  the  region  does  not  approach  zero 
but  instead  fills  in,  developing  a  shape  of  its  own. 


Figure  10.  Field  for  square  transducer  (lossy  case) . 

Another  major  difference  is  discovered  when  comparing 
Figures  9  and  11.  Where  the  initial  amplitude  for  the  loss¬ 
less  case  was  1.0,  the  amplitude  of  the  spatial  excitation, 
the  initial  amplitude  for  the  lossy  case  has  been  reduced  to 
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approximately  0.30.  This  reduction  in  amplitude  was  ob¬ 
served  for  all  test  cases  with  a  linear  loss  coefficient. 


The  results  obtained  using  a  circular  transducer  are 
similar  to  those  obtained  with  the  square  transducer.  The 
spatial  excitation  for  a  circular  transducer  of  width  15  and 
amplitude  1.0  is  shown  in  Figure  12.  Here  a  width  of  15  was 
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Figure  11.  Side  view  for  lossy  case. 


used  to  improve  the  resolution  of  the  circular  excitation. 
The  spatial  impulse  response  using  this  source  is  shown  in 
Figure  13.  Comparison  of  Figures  8  and  13  illustrates  the 


differences  between  the  spatial  excitation  shapes.  The  ma¬ 
jor  difference  is  the  tails,  which  now  have  a  circular 
cross-sectional  shape.  Some  minor  variations  can  also  be 
observed  in  the  large  initial  response  early  in  time. 
Figure  14  shows  the  spatial  impulse  response,  using  the 
circular  excitation  of  Figure  12,  for  the  lossy  case.  The 
previous  comments  for  the  square  transducer  concerning  the 


Figure  12.  Circular  piston  spatial  excitation. 

difference  in  amplitudes  and  shape  of  the  "tails''  for  the 
lossless  and  lossy  cases  also  apply  to  the  circular 


transducer.  Figure  15  is  a  side  view  plot  of  Figure  14 
showing  the  reduction  in  amplitude  is  also  present  with  the 
circular  transducer. 

The  results  obtained  when  using  the  square-piston  and 
circular-piston  spatial  excitation  for  both  lossless  and 
lossy  media  compare  favorably  with  those  obtained  by  Guyomar 
and  Powers  [Refs.  1-3,9].  Based  on  these  results,  it  was 
determined  that  the  program  correctly  implemented  the 
propagation  models,  and  therefore  could  be  used  to  compute 
the  spatial  impulse  response  using  new  types  of  excitation. 


B.  NEW  EXCITATION  SHAPES 

The  program  has  to  ability  to  analyze  any  excitation 
source  shape  that  will  fit  on  the  64  x  64  base  array.  Two 
new  single  excitation  shapes  and  four  multiple  element  exci¬ 
tation  shapes  are  in  the  program .  The  results  for  three  of 


Figure  14.  Field  for  circular  transducer  (lossy  case) . 

these,  the  pyramid,  the  5  element  linear  array  with  constant 
amplitude,  and  the  5  element  linear  array  with  stepped  am¬ 
plitude  are  included. 
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A  pyramid  shaped  spatial  excitation  of  width  11  and  am¬ 
plitude  1.0  is  shown  in  Figure  16.  The  computed  spatial  im¬ 
pulse  response  for  this  excitation  in  lossless  media  is 
shown  in  Figure  17.  Since  this  particular  waveshape  has  a 
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Figure  15.  Side  view  for  lossy  case. 


low  spatial  frequency  content,  the  shape  of  the  wave  stays 
much  the  same  for  small  values  of  time.  Figure  18  shows  the 
spatial  impulse  response  for  the  same  pyramid  excitation  in 


Figure  19  shows  the  5-element  linear  array  configura¬ 
tion.  For  this  array  the  elements  are  each  5x5  with  their 
centers  positioned  10  units  apart.  This  produces  a  space  of 
5  units  between  each  element.  The  calculated  spatial  im¬ 
pulse  response  for  this  excitation  is  shown  in  Figure  20. 


Figure  16.  Pyramid  shaped  spatial  excitation. 

The  shape  for  each  individual  excitation  is  similar  to  the 
single  excitation  response  in  Figure  8  until  the  tails  start 
to  spread  out  and  interfere  with  each  other.  In  the  lossy 
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case.  Figure  21,  the  pattern  is  similar  but  with  increased 
interference  in  the  region  between  the  individual  responses. 


Figure  17.  Field  for  pyramid  excitation  (lossless  case). 

The  spatial  excitation  for  the  last  test  case  presented 
is  shown  in  Figure  22.  This  is  a  5-element  linear  array 
similar  to  Figure  19  but  with  stepped  amplitude.  The  center 
element  has  amplitude  of  1.0  as  before.  The  two  elements 
adjacent  to  the  center  have  amplitudes  of  2/3,  and  the  outer 
two  elements  have  amplitudes  of  1/3.  This  particular  case 
is  presented  to  show  the  flexibility  of  the  program.  Figure 


23  shows  the  spatial  impulse  response  in  lossless  media. 
The  results  are  similar  to  those  of  Figure  20  except  for  the 
dominance  of  the  center  element.  Figure  24  shows  the  spa¬ 
tial  impulse  response  for  the  same  excitation  in  lossy  me¬ 
dia. 


Figure  20.  Field  for  five  element  linear  array 

(lossless  case) . 


Figure  23.  Field  of  five  element  linear  array  with 
stepped  amplitude  (lossless  case) . 
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Figure  24.  Field  for  five  element  linear  array  with 
stepped  amplitude  (lossy  case) . 
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V.  SUMMARY 


This  thesis  has  investigated  the  ability  to  perform 
complex  simulations  of  scalar  wave  propagation  through  loss¬ 
less  media  and  media  with  a  loss  coefficient  linear  in  fre¬ 
quency  on  a  micro-computer.  A  program  was  written  in 
Microsoft  Fortran  V3.31  which  allows  the  user  to  select  or 
specify  the  initial  conditions  and  media  type,  then  calcu¬ 
lates  the  resulting  spatial  impulse  response.  The  set  goal 
of  thirty  minutes  for  each  simulation  was  achieved  through 
code  optimization  and  partial  array  calculation  permitted  by 
symmetry . 

After  initial  program  verification,  the  spatial  impulse 
responses  for  both  square-piston  and  circular-piston  spatial 
excitation  were  computed.  The  results  obtained  agreed  with 
previously  accepted  results.  The  spatial  impulse  response 
was  then  calculated  for  three  new  spatial  excitation 
sources,  a  pyramid-shaped  excitation,  a  five  element  linear 
array  excitation  with  constant  amplitude,  and  a  five  element 
linear  array  excitation  with  stepped  amplitude.  The 
computed  spatial  impulse  response  for  all  simulations  have 
been  plotted  and  are  presented  for  visual  interpretation. 

Future  development  should  concentrate  in  two  areas. 
First,  the  third  model  for  media  with  a  loss  coefficient 
which  is  quadratic  in  frequency  must  be  added  to  the 


program.  This  will  provide  all  the  tools  necessary  for 
analysis  of  wave  propagation  through  air,  liquid,  gas,  and 
biological  tissue.  Second,  the  base  array  size  must  be  in¬ 
creased  beyond  the  present  64  x  64.  An  array  size  of  64  x 
64  proved  adequate  for  program  development  and  testing  but 
increased  size,  and  therefore  increased  resolution,  is  nec¬ 
essary  for  serious  analysis  and  evaluation.  The  limitation 
in  all  presently  available  micro-computers  which  restricts 
array  size  is  processing  speed.  For  example,  increasing  ar¬ 
ray  size  to  128  x  128  would  increase  execution  time  by  a 
factor  of  four,  requiring  almost  two  hours  for  each  simula¬ 
tion  run.  Transferring  program  development  to  one  of  the 
newer  micro-computers,  utilizing  the  Intel  80386  micro¬ 
processor  will  allow  larger  array  sizes  while  maintaining 
reasonable  execution  times. 


This  appendix  contains  the  FORTRAN  source  code  for  the 
program  developed  in  the  thesis.  This  program  was  compiled 
and  run  using  Microsoft  Fortran  Version  3.31. 


3LARGE 

SNOFLOATCALLS 

C 

c 

C  COMMENTS: 

C 

C  THIS  PROGRAM  HAS  WRITTEN  BY  LT.  TIMOTHY  *RRILL.  USB 

C 

C  IT  is  WRITTEN  IB  MICROSOFT  FORTRAN  73. 1  i  AND  IS  CAPABLE  OF  RUBBING 

C  OB  ANY  IBM  OR  COMPATIBLE  MICRO -COMPUTER  WITH  LITTLE  OR  BO 

C  MODIFICATION .  THE  USE  OF  A  8087  NUMERICAL  CO-PROCESSCR  IS  HIGHLY 
C  RECOFMENDED  TO  REDUCE  EXECUTION  TIME. 

C 

C  TO  REDUCE  EXECUTION  TIME  THE  PROGRAM  IS  WRITTEN  TO  USE  A  8*  X  8* 

C  BASE  ARRAY  SIZE.  AND  IS  WRITTEN  AS  SINGLE  PRECISISOR.  FUTURE 
C  MODIFICATION  TO  A  128  X  128  (GR  LARGER)  ARRAY  SIZE  AND/OR  DOUBLE 

C  PRECISION  IS  VERY  EASY. 

C 

C  THE  PROGRAM  CALLS  THE  FOLLOWING  SUBROUTINES  WHOSE  SOURCE  CODE  IS 
C  NOT  INCLUDED  HERE: 


C  FFT2C  IMSL  ROUTINE  TO  COMPUTE  FFT  OF  A  COMPLEX 

C  VALUED  SEQUENCE  OF  LENGTH  EQUAL  TO  POWER 

C  OF  TWO 

C  FFT3D  IMSL  ROUTINE  TO  COMPUTE  THE  FFT  OF  A  COMPLEX 

C  VALUED  1.2,  OR  3  DIMENSIONAL  ARRAY 

C  PFTCC  IMSL  ROUTINE  CALLED  BY  FFT3D 

C  MMBSJN  IMSL  ROUTINE  TO  CALCULATE  BESSEL  FUNCTION  OF 

C  THE  FIRST  KIND  OF  ZERO  ORDER  WITH  REAL  ARGUMENT 

C  MMBSIN  IMSL  ROUTINE  TO  CALCULATE  MODIFIED  BESSEL 

C  FUNCTION  OF  THE  FIRST  KIND  OF  ZERO  ORDER  WITH 

C  REAL  ARGUMENT 

C 

C  T.D.M.  20  JANUARY  1987 

C  *»*»**»»«*******«***.****»***..»***»**...***»*...»****».*.»**»**.* 

c 

REAL  A. A2.B.B1, BASE. BI.BR.C.C2, CENTER 


c 


B«AT 

MAI. 

R EAL 
REAL 

COMPLEX 
COMPLEX 
C 

INTEGER 

urara 

integer 

c 

DIMENSION 
DIMENSION 
DIMENSION 
C 

CHARACTER 
CHARACTER 
CHARACTER 

C 

C  ****«**»«, 

C  CONSTANTS 

£  ********************************************* A** 4*******t******** 

SPEED  OF  SOUND  IN  VACUUM  (1500  METERS/SEC) 

SPEED  OF  SOUND  SQUARED  (METERS**2/SEC**2) 

HALF  OF  NROWS  (OR  NCOLS)  PLUS  1 
INTEGER  USED  TO  SPECIFY  SIZE  OF  VECTOR  SENT  TO 
FFT2C  (REFER  TO  IMSL  SOURCE  CODE  FOR  FFT2C) 

USED  FOR  GAUSSIAN  IMPULSE  FUNCTION 
USED  ON  BESSEL  FUNCTION  CALLS  TO  SIGNIFY  ZERO  ORDER 
DEFAULT  NUMBER  OF  COLUMNS 
DEFAULT  NUMBER  OF  ROWS 
USED  FOR  GAUSSIAN  IMPULSE  FUNCTION 
SIGMA  SQUARED 

NUMBER  OF  RCMS  SET  TO  ZERO  -  SIMULATES  STEP  FUNCTION 
2  *  PI 

C  A******************************************.**********.******.... 
C  -  1500.0 
N  -  1 
M  -  0 
STEP  -  4 
NROWS  -  64 
NCOLS  -  64 
PI  -  3.14159265 
SIGMA  -  1.0 
MEANX  -  0.0 
C2  -  C  *  C 
H  -  NROWS/ 2  +  1 
TWOPI  -  2.0  *  PI 
SIGMA2  -  SIGMA  *  SIGMA 
C 


C  C 

C  C2 

C  H 

C  M 

C 

C  MEANX 

C  N 

C  NCOLS 

C  NROWS 

C  SIGMA 

C  SIGMA2 

C  STEP 

C  TWOPI 


BOOTH ,  IRC ,  MAXTIME ,  MANX ,  PI 

R123,  RESULT, RB01.RB02.RMK,  SIGMA.  SIGMA2 

TEMP1.TZMP2,  TEMPS,  TEMP4.TEMP3 

TIME , IS ,  TWOPI , X , Y , Z , Z2 , ZPRIME , TVECTR 

C123 .CWK.CTEMPl .CTEMP2 .CTEMP3 
PULSE, WORK 

COL , FTYPE , H , I , IER , I JOB , IWK , IWK1 , J 

M.H.  HAtCEN,  NCOLS.  NROWS,  R1.R2 

ROW , RBO , SHAPE , STEP , WIDTH , XPOS , YPOS 

CWK(64 ) , IWK( 534 ) . IWK1 ( 7 ) , RHOl ( 64 ) , RWK< 534 ) 
PULSE ( 64 . 64 ) , RESULT ( 64 , 84 ) , RB02 ( 6 4 , 64 ) 

TIME ( 64 ),  WORK ( 64 , 64 ), ZPRIME ( 64 ) 

FN»8,FHAME*8,T1,T2 

EXT1*4 , EXT2*4 , EXT3*4 , EXT4*4 

FNAME1*12,  FNAME2*12,  FNAME3*12,  FNAME4*12 
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3ST  USER  INPUT 

EPOS  -  1 
TEOS  -  2 
CALL  CLRSCR 
CALL  OOTCQCYd,  1) 

‘BITER  REQUESTED  VALUES.  RECOMMENDED  VALUES  IN  (  )• 
CALL  OOTCKTCTPOS.XPOS) 

WRITE(*,*)  'ENTER  FILENAME  (8  CHAR  HZ'?):  * 

CALL  GOiaXY(YPOS,X?OS+32) 

READ(*, ’ (A)* )  PN 
CALL  GOTOXY(YPOS+l,XK>S) 

HRITK(*.*)  ‘SELECT  PILTER  TYPE’ 

CALL  GOTOXY(YPOS+2,XFOS+8) 

WRITE(* ,*)  ><1»  GP1  (LOSSLESS  MEDIA)1 
CALL  OOTC3XY  ( YFOS+3 ,  XPOS-t-8 ) 

WRITE(*,»)  ’ <2>  OP2  (LOSSY  MEDIA) 1 
CALL  GOTOXY(YPOS+4 ,XPQS) 

WRITE(*,9)  ■  ENTER  1  OR  2:  ' 

READ(*,»)  PTTPE 

CALL  GOTOKY(YPOS+3,XFOS) 

VEUTEC*,*)  ‘SELECT  SOURCE  SHAPE' 

CALL  OOTCKY(YPOS+C,XPOS+* ) 

VRITEC*,*)  ‘«1>  SQUARE' 

CALL  QOTOXY( YPOS+7 , XPOS+4 ) 

WUTE(*.*)  ‘<2>  CIRCLE* 

CALL  GOTOXY(YPOS+8,XFOS+*) 

WUTE(«,*)  ‘«3»  GAUSSIAN’ 

CALL  QOTOXY( YFOS+9 .XPOS+4 ) 

WUTE(*,*)  '«♦>  PYRAMID' 

CALL  GOTCDCYt  YPOS+10  .XPOS+4  ) 

WRITEC*.*)  ‘<5>  TRUNCATED  PYRAMID’ 

CALL  aOTOKY(YPOS+ll.XPOS+A) 

V*UTE(«,«)  *<8»  3  PULSE  LINEAR  ARRAY  (SQ.)‘ 

CALL  OOTOXY(YPOS+1.2,XPOS+*) 

WRITE(*,*)  ’<7»  3  PULSE  LINEAR  ARRAY  (GAUSS.)’ 

CALL  GOT0KY(YPOS+8.XPO3+*0) 

WRITE(*.*>  '  <8>  5  PULSE  LINEAR  ARRAY  (STEPPED)’ 

CALL  GOTOXY(YPOS+7.XPOS-H»0) 

V®ITE(*.*)  ’  «9>  3  PULSE  LINEAR  ARRAY  (PARABOLA)’ 

CALL  OOTOXY( YPOS+8 , XPOS+40 ) 

WRITE(«,*)  ><10>  RESERVED  FOR  FUTURE  USE’ 

CALL  GOTOKY( YPOS+9 , XPOS+40 ) 

WUTE(*,*)  ’«U>  RESERVED  FOR  FUTURE  USE’ 

CALL  GOTOXY(YPOS+10,XPOS+*0) 

VNUTE(*,»)  ’<12>  RESERVED  FOR  FUTURE  USE’ 

CALL  GOTOXY(YPOS+11,XPOS+40) 

HRITE(*,»)  ’<13>  RESERVED  FOR  FUTURE  USE' 

CALL  aOTOKY(YFOS+12,XFOS+40 ) 

WUTE(*,*)  ’<!*»  RESERVED  FOR  FUTURE  USE’ 
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cm  <xmacT(TK*+u,xpo6) 

WZTI<*.8>  '  ENTER  SELECTION: 

READ(*,«)  SHAPE 
IF  (SHAPE  .LT.  6)  THEM 

Cm  OOTCOCY(YP06+1*,XPOS) 

WRITE(*,9)  '  EMTER  SOURCE  SIZE  (ODD  INTEGER  <-  63)  (11): 
READ(*,*)  WIDTH 

EMDIP 

Cm  GOT0KY{YFO8+15,XFO8) 

MUTE(*.*)  ’ENTER  RHO  (200):  ’ 

Cm  GOTOXY(Y?OS+lJ,XH)S+2S) 

READ(*.*)  RBO 
Cm  aOI«Y(YPOS+16,XPOS) 

WITE(*,*)  ’ENTER  Z  (METERS)  (0.1):  ’ 
cm  00Tacr(TP0s+i6,xpcs+26) 

RIAD(*,*)  Z 

CALL  QOTCKY  ( YPOS+ 1 7 ,  XPCS ) 

HRITE(*,7)  Z/C 

CALL  GOTOXY(YPOS+17 .XPOS+M ) 

READ<*,»)  MAXTIME 
IF  (FTTFE  .EQ.  2)  THEN 

CALL  GOTOXYCYPOS+IS.XFOS) 

MUTE(*,*)  ’ENTER  TALUK  FOR  A  (0.008):  ’ 

Cm  GOTaXY(YFOS+18,XFOS+28) 

READ(*,*)  A 
A2  -  A  *  A 
ESDI? 

FORMAT ( ’  ENTER  MAXTIME  (REAL  >  ’ ,  1PE1A.6*,  ’  ) :  *  ) 

FORMAT  (A\ ) 

HWIDTH  -  REAL  (WIDTH)./ 2 
Z2  -  Z  •  Z 
TS  -  Z/C 
cm  CLRHCR 
cm  aoTOKTd.i) 

HRITE(*,*)  ’CALCULATING...’ 

INITIALIZE  ARRAYS  AND  VECTORS  TO  ZERO 


DO  3  ROW  -  1, KNOWS 
TDC(ROW)  -  0.0 
RHOl(ROW)  -  0.0 
DO  3  COL  -  1.NC0LS 

PULSE (COL, ROW)  -  (0. 0.0.0) 
WORX(COL.ROW)  -  (0.0. 0.0) 
RESULT (COL, ROW)  -  0.0 
R8D2(C0L,R0W)  -  0.0 
CONTINUE 


OPEN  OUTPUT  FILES 


!>»«««« 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 


PROGRAM  SETS  SAME  OF  PROBLEM  <FN>  FROM  USER  AND  FORMS: 
FILENAME. EXT  FILE  DESCRIPTION 


<FR>.  IMP  ARRAY  OF  INPUT  IMPULSE  FUNCTION 
<FN».TXT  SWMARY  INFORMATION  FOR  THIS  PROBLEM 

<FN>.DAT  FINAL  RESULTS  ARRAY 

ROTE  1:  ALL  FILES  ARE  IN  ASCII  TEXT  FORMAT 

EXT1  -  '.IMP' 

EXT2  -  ' .TXT' 

EXT3  -  ’ .DAT' 

RAMLEN  -  0 
T1  -  '  ’ 

00  10  1-1,8 
T2  -  FN(I:I) 

IF  (T2  .HE.  Tl)  THEN 
FHAME(I.-I)  -  FN(I : I ) 

HAMLEN  -  NAMLEN  +  1 
END  IF 
CONTINUE 

FNAMEK1:  NAMLEN)  -  FRAME 
FRAME  KHAMLES+l :  HAMLEN+A )  -  EXT1 
FRAME2(1: NAMLEN)  «  FRAME 
FHAME2 ( RAMLER+1 : NAMLER+* )  -  EXT2 
FRAME3 ( 1 : RAMLER )  -  FRAME 
FNAME3(NAtUR+l:RAMLEN+A)  -  EXT3 
IF  (HAMLEN  .LT.  8)  THEN 
FHAMEKRAMLER+5: 12)  -  ' 

FRAME2(NAMLEN+3: 12)  -  ’ 

FRAME3 (RAMLER+3 : 12 )  -  ’ 

END  IF 

OPEN ( 1 , FILE-FHAME1 . STATUS- ' NEW ' ) 

OPEN ( 2 , FILE— FNAME2 , STATUS— ' HEW ' ) 

OPEN  (  3 ,  FILE-FNAME3 ,  STATUS-  *  HEW ' ) 


WRITE  HEADER  INFORMATION  TO  «FN>.TXT 


WRITE (2, 11)  FR 
IF  (FTYPE  ,EQ.  1)  THEN 

WRITE(2,*)  ’  FILTER  IS  GP1  FOR  LOSSLESS  MEDIA' 

ELSEIF  (FTYPE  .EQ.  2)  THEN 

WRITE(2, *)  '  FILTER  IS  GP2  FOR  LOSSY  (LINEAR)  MEDIA' 
ENDIF 

IF  (SHAPE  .EQ.  1)  THEN 

WRITE (2, 12)  WIDTH, NCOLS.NROWS 
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ELSEIF  (8BA7X  .IQ.  2)  TIBI 

MUTXC2, 13)  WIDTH, NCOLS,*ROWB 
ELSEIF  (SHAW  .IQ.  3)  THU 

MRITE(2.14)  WIDTH, WOOLS, BROWS 
ELSEIF  ! SHAPE  .IQ.  *)  THU 

tStITE(2, 111)  WIDTH ,  HCOLS ,  HRONS 
ELSEIF  (SHAR  .IQ.  3)  THE* 

WHITE (2. 112)  WIDTH, NCOLS, HROWS 
ELSZIT  (SHAR  .IQ.  6)  THE* 

(CITE(2,113) 

ELSEIF  (SHAR  .IQ.  7)  THE* 

HRITK2, 114) 

ELSEIF  (SHAR  .EQ.  8)  THE* 

VSITE(2,113) 

ELSEIF  (SHAR  .EQ.  9)  THE* 

HRITE(2,116) 

ERDIF 

WRITE (2, 15)  STEP 
WRITE (2, 16)  Z 
WRITE (2, 17)  RHO 
WRITEC2, 18)  MAXTIME 
IF  (FTTR  .EQ.  2)  THEN 

WRITE(2, 19)  A  ’ 

EKDIF 

FORMAT! 'SUPMARY  INFORMATION  FOR  ’  ,A /) 

FORMAT ( ’  IMPULSE  FUNCTION  IS  A  SQUARE  WITH  WIDTH  ',12, ’  ON  A  HAS 

♦E  OF', 13,*  X* ,13) 

FORMAT (’  IMPULSE  FUNCTION  IS  A  CIRCLE  WITH  DIAMETER  ’,12, '  ON  A 

♦BASE  OF',13,'  X'. 13) 

FORMAT ( '  IMPULSE  FUNCTION  IS  GAUSSIAN  WITH  DIAMETER  ',12,'  ON  A 

♦EASE  OF',13,*  X* ,13) 

FORMAT ( '  IMPULSE  FUNCTION  IS  PYRAMID  WITH  WIDTH  ’,12,'  ON  A  BAS 
+E  OF',13,'  X' .13) 

FORMAT ( ’  IMPULSE  FUNCTION  IS  TRUNCATED  PYRAMID  WITH  WIDTH  ',12,' 
+ON  A  BASE  OF',13,’  X',I3) 

FORMAT ('  IMPULSE  FUNCTION  IS  5  PULSE  LINEAR  ARRAY  WITH  SQUARE  PU 
+LSES') 

FORMAT ('  IMPULSE  FUNCTION  IS  3  PULSE  LINEAR  ARRAY  WITH  GAUSSIAN 
♦PULSES') 

FORMAT! '  IMPULSE  FUNCTION  IS  3  PULSE  LINEAR  ARRAY  WITH  STEPPED  A 
+WPLITUDE' ) 

FORMAT! '  IMPULSE  FUNCTION  IS  5  PULSE  LINEAR  ARRAY  WITH  PARABOLIC 
+  AMPLITUDE') 

FORMAT! '  STEPS  I ZE  -  ’,12) 

FORMAT! '  Z  -  ’,1PE14.7,’  METERS') 

FORMAT ( ’  RHO  -  ' ,13) 

FORMAT!  ’  MAXTIME  -  ',1PE1<..7,'  SECONDS’) 

FORMAT!'  A  -  ’ .1PE14.7) 


rl  N  v 4  A. 


c 

c 

c 


INITIALIZE  IM  VECTOR 


20 

C 

C 

C 

C 


21 

22 

C 

C 

c 

c 

c 

c 


30 

c 

c 

c 

c 

c 

c 


30 

c 


60 

C 


70 

C 


INC  -  (MAXTIME-TS  )  /REAL  (NRQWS-STEP- 1 ) 
00  20  I  -  <STEP+l),NROWB 

TOK(I)  -  TS  +  (I-STZP-1)  *  INC 
CORT I  VUE 


OUTPUT  TIME  SUH4ARY  TO  <FN>.TXT 


HRITE(2,*) 

VSXTE(2,9>  ’TIME  SUMMAXY:' 
t*ITE<2.21)  T8.STEP+1 
WRITE (2, 22)  IRC 

FORMAT ( ’  ZERO  TIME  STARTS  AT  T-Z/C  —.1PE14.7,’  SECONDS  IK  RQWC 

+.12,’)*) 

FORMAT ( ’  AND  INCREMENTS  BY  \1PEU.7;’  SECONDS  PER  ROW’/) 


CREATE  Z-PRIME  VECTOR:  SQUARE  EACH  VALUE  OF  TIME, 

MULTIPLY  BY  C  SQUARED  AND  SUBTRACT  Z  SQUARED.  THEN 
TAKE  SQUARE  ROOT. 

******w***fc*******************«********«************************* 
DO  30  I  -  ( STEP+2 ) , NROW3 

ZPRIME(I)  -  SQRT( (TIME(I )  *  TIME(I)  •  C2)  -  Z2) 

CONTINUE 


•••♦•*••*******•*•*»**•****•••*•****»*»•*»*»*******•*********»**» 
INITIALIZE  RHO  VECTOR  AND  GENERATE  33  X  64  RHO  MATRIX 
ONLY  ONE  HALF  OF  THE  MATRIX  HAS  TO  BE  CALCULATED  SINCE 
THE  OTHER  SIDE  IS  SYMMETRICAL 

***************************************************************** 
INC  “  REAL (RBO) /REAL (NCOLS/ 2-1) 

BASE  -  REAL(-RHO)  -  INC 
RBOl(l)  -  BASE 
DO  30  I  -  1 , NCOLS- 1 

RBOld+1)  -  BASE  +  I  *  INC 
CONTINUE 

DO  60  I  -  l.NROWS 

RB02CH.I)  -  RHOld) 

RH02(I,H)  -  RHOld) 

CONTINUE 

DO  70  ROW  -  l.NROWS 
DO  70  COL  -  1,B 

RH02(C0L,R0W)  -  SQRT(RH02(C0L,H)«*2  +  RH02(H ,ROW)»*2) 
CONTINUE 
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c 

c 


OUTPUT  RBO  SUMAXY  10  «FN».TXT 


41 

42 

43 

44 


WIR(2,0)  'SBO  SUPMARY: ' 

MUR(2.44)  H.HCOLS/2+1 
WITXU.Al)  BASE 
WITS  (2, 42)  ISC 
MUTE(2,43)  RB02(1,1) 

FORMAT!  ’  BASE  VALUE  -  (-RBO  -  ISC)  -  MPE14.7) 

FORMAT ( ’  INCREMENT  -  RH0/(64/2  -  1)  -  \1PE14.7) 

FORMAT (  ’  MAXIMUM  VALUE  IS  AT  RB02U.1)  -  \XPE14.7) 

FORMAT! '  ARRAY  IS  CENTERED  ABOUT  (‘,13,',’, 13 ,')') 


C 

C 

C 

C 

C 

C 


1200 

C  ‘ 

C 


1300 

C 

C 


1400 

C 

C 


GENERATE  EXCITATION  -  CENTER  POINT  MUST  BE  AT  (33,33) 


«*  GENERATE  SQUARE  EXCITATION  *** 

IF  (SHAPE  .EQ.  1)  THEN 

DO  1200  ROW  -  (NROWS/2+1-WIDTH/2), (NR0WS/2+1+WIDTH/2) 

DO  1200  COL  -  (NCOLS/2+1-WIDTH/2) , (NC0LS/2+1+WIDTH/2) 

PULSE (COL, ROW)  -  (1.0, 0.0) 

CONTINUE 

«•  GENERATE  CIRCULAR  PULSE  *** 

ELSEIF  (SHAPE  .EQ.  2)  THEN 

DO  1300  ROW  -  (NROWS/2+1-WIDTH/2), (NROWS/2+1+WIDTH/2) 

DO  1300  COL  -  (NCOLS/2+1-WIDTH/2) , (NC0LS/2+14WIDTH/2) 

IF  (9QRT(REAL(COL-H)**2+REAL(ROW-H)**2)  .LT.  HWIDTH)  THEN 
PULSE (COL. ROW)  -  (1. 0,0.0) 

END  IF 
CONTINUE 

**»  GENERATE  CIRCULAR  EXCITATION  WITH  GAUSSIAN  AMPLITUDE  *** 
ELSEIF  (SHAPE  .EQ.  3)  THEN 

TEMPI  -  l/SQRT(TWOPI*SIGMA2) 

TEMP2  -  1/ (2*SIGMA2) 

DO  1400  ROW  -  (NROWS/ 2+1- WIDTH/2) , (NROWS/2+1+WIDTH/2) 

DO  1400  COL  -  (NCOLS/2+1-WIDTH/2) , (NCOLS/2+1+WIDTH/2) 

TEMP3  -  SQRT(REAL(COL-H)**Z+REAL(ROW-H>«*Z) 

IF  (TEMPS  .LT.  HWIDTH)  THEN 

PULSE  (COL,  ROW)  -  TBMP1*EXP(-TEMP2*(TEMP3-MEANX)**Z) 

END  IF 
CONTINUE 

•**  GENERATE  PYRAMID  EXCITATION  *** 

ELSEIF  (SHAPE  .EQ.  4)  THEN 
R1  -  (NROWS/2)  -  (WIDTH/2) 

R2  -  (NROWS/2  +  2)  +  (WIDTH/2) 

INC  -  1.0 /(WIDTH/ 2  +  1) 
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mass* 


<  v.'/-  X  •  "  •/  1 

/JaWaaV. 


DO  1900  I  -  1, WIDTH/2 

DO  1900  BOW  -  R1+I.H2-I 
DO  1900  COL  -  R1+I.R2-I 
PULSE (COL, BOM)  «  INC  *  I 

1900  CONTINUE 

PULSE (H,H)  -  1.0 
C 

C  **»  GENERATE  TRUNCATED  PYRAMID  EXCITATION  **• 

ELSSIF  (SHAPE  .EQ.  9)  THEN 
R1  -  (NBOHS/2)  -  (WIDTH/2) 

R2  -  (NBOWS/2  +  2)  ♦  (WIDTH/2) 

INC  -  1.0/ (NCOLS/2) 

BASE  -  (NCOLS/2  -  WIDTH/2  +  1)  *  INC 

INC  -  1.0 /H 

DO  1800  I  -  1, WIDTH/2 

DO  1600  BOW  -  R1+I.R2-I 
DO  1600  COL  -  R1+I.R2-I 

PULSE (COL. ROW)  -  BASE  +  INC  *  I 

1600  CONTINUE 

PULSE(H.B)  -  1.0 
C 

C  ***  GENERATE  S  ELEMENT  LINEAR  ARRAY  -  SQUARE  *** 

ELSEIF  (SHAPE  .EQ.  S)  THEN 
WIDTH  -  9 

DO  1700  I  -  1, WIDTH 
COL  -  I  *  10 
DO  1700  J  -  1, WIDTH 
COL  -  COL  +  1 

DO  1700  ROW  -  (H  -  WIDTH/2), (H  +  WIDTH/2) 

PULSE (COL, ROW)  -  (1.0, 0.0) 

1700  CONTINUE 

C 

C  *»*  GENERATE  5  ELEMENT  LINEAR  ARRAY  -  GAUSSIAN  *** 
ELSEIF  (SHAPE  .EQ.  7)  THEN 
WIDTH  -  5 

TEMPI  -  1 / SORT ( TW0PI*SIGMA2 ) 

TEMP2  -  1/ (2*SIGMA2) 

DO  1800  ROW  -  (NR0WS/2+1-WIDTH/2) , (NRCWS/2+1+WIDTH/2) 

DO  1800  COL  -  (NCOLS /2+1-WIDTH/ 2) , (NC0LS/2+1+WIDTH/2) 
TEMP3  -  SQRT (REAL (COL-H)**2+REAL(ROW-H)**2) 

IF  (TEMP  .LT.  REAL(WIDTH/2) )  THEN 

TEMP*  -  TEMP1*EXP(-TEMP2*(TEMP3-MEANX)**2) 
PULSE(COL-20 .ROW)  -  TEMP* 

PULSE (COL- 10, RCN)  -  TEMP* 

PULSE ( COL, RCM)  -  TEMP* 

PULSE (COL+IO, ROW)  -  TEMP* 

PULSE (COL+20, ROW)  -  TEMP* 

ENDIF 

ia  oo  ccmimt 

c 
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c  ***  generate  5  zlbcvt  linear  array  -  stepped  *** 


WIDTH  -  9 

n»i  -  1. 0/3.0 

TWF2  -  2. 0/3.0 

TEMP3  -  1.0 

DO  1000  COL  -  1,  WIDTH 

DO  1000  BOW  -  (H-WIDTH/2) . (H+WIDTH/2) 

PULSE (COL+IO, ROW)  -  TOffl 
PULSE (COL+20, BOW)  «  TEMP2 
PULSE (COL+30. ROW)  -  TEMP3 
PULSE (COL+*0, ROW)  -  TEMP2 
PULSE (COL+50 .ROW)  -  TSffl 
1000  continue 
c 

C  ***  LINEAR  ARRAY  OP  3  ELEMEHTS  -  PARABOLIC  AMPLITUDE  *** 

C  IMPLEMENTS  THE  FOLLOWING  FORMULA: 

C  (X  -  H)**2  -  -*  *  P  *  (Y  -  X) 

C  WITH 

C  H  -  0 

C  P  -  236 

C  K  -  236 

C  TO  GIVE  A  DOWNWARD  OPENING  PARABOLA  THAT  HAS  A  PEAK  OF  1.0 

C  AMD  IS  EQUAL  TO  ZERO  AT  X  -  -32.  32  (ROWS  2  AND  6*). 

C  EQUATION  IS  EVALUATED  AT  X  -  0,  10  AND  20. 

C 

ELSE IF  (SHAPE  .EQ.  0)  THEN 
WIDTH  -  5 

TEMPI  -  (*.*236.  -  20**2 ) / ( * . *256 . ) 

TEMP2  -  (*.*236.  -  10**2)/ (*.*236. ) 

TEMP3  -  (*.*236.  -  0**2)/(*.*236. ) 

DO  2000  COL  -  1, WIDTH 

DO  2000  ROW  -  (H-WIDTH/2). (H+WIDTH/2) 

PUL3E(COL+10,ROW)  -  TEMPI 
PULSE (COL+20, ROW)  -  TEMP2 
PULSE(COL+30,RQW)  -  TEMP3 
PULSE (COL+*0, ROW)  -  TEMP2 
PULSE (COL+30, ROW)  -  TEMPI 
2000  CONTINUE 
C 

ENDIF 

C 

c  *******«.*••*•*****••*••*.***•**»***•••*****•*•#•*•***•***»**• 

C  OUTPUT  IMPULSE  FUNCTION  ARRAY  TO  <FN>.IMt 

DO  3000  I  -  l.NROWS 

WRITZ(  1 , 102)  (REAL(PULSE(J.D).J  -  l.NCOLS) 

3000  CONTINUE 
CLOSE! 1) 

C 
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PASS  PULSE  TO  2-D  FFT  SUBROUTINE.  I  JOB  -  1  (CARS  TAKE  FFT 


I  JOB  -  1 

CALL  FFT3  D(  PULSE  ,NCOL3,  NROWS , NCOLS ,NRCW8 , 1 , 1 JOB , IWK.RWK.CWK) 


SHIFT  RESULTS  -  IE:  SWITCH  QUADRANTS  1  AND  3,  2  AND  4  IK  THE 
COL, ROW  PLAHE 

A************************************************************* i 

CALL  SHIFT  (PULSE, RROWS.MCOLS) 


HAIR  PROGRAM  STARTS  HERE 

CALCULATIONS  START  AT  ROW  -  STB P+1. 

ROTE  THAT  ORLY  ORE  QUARTER  OF  THE  WORK  MATRIX  HAS  TO 
BE  CALCULATED.  THE  REMAINING  THREE  QUARTERS  IS  FORMED  BY 
FOLDING  AND  FLIPPING.  THIS  IS  POSSIBLE  DUE  TO  SYMMETRY. 
************************************************************ , 

CALL  GOTQXYd,  1) 

WRITE(»,»)  -CALCULATING  ROW  ' 

IP  (FTYPE  .EQ.  1)  THEN 
DO  200  I  -  STEP+1 , NROWS 
CALL  GOTQXYd, 18) 

WRITE<*,101)  I 

DO  210  ROW  -  1,H 
DO  210  COL  -  1,B 

B1  -  ZPRIME(I)  *  RH02(C0L.R0W) 

CALL  (MBSJHCBl.N, TEMPI, 1ER) 

WORK (COL, ROW)  -  TEMPI  *  PULSE ( COL , ROW) 

‘CONTINUE 


FLIP  QUADRANT  2  TO  QUADRANT  3 


J  -  0 

DO  220  ROW  -  H+l, NROWS 
J  -  J  +  2 
DO  220  COL  -  l.H 

WORK (COL , ROW)  -  WORK(COL.ROW-J) 
CONTINUE 


FLIP  LEFT  SIDE  TO  RIGHT  SIDE 


J  -  0 

DO  213  COL  -  H+l.NCOLS 
J  -  J  +  2 

DO  215  ROW  -  1, NROWS 

WORX(COL.ROW)  -  WORK(COL-J.ROW) 
CONTINUE 


I,  .+»».  w. 


j»  -<fc  ,<>  J>  »'• 


TAKE  DIVERSE  ITT  OT  EACH  COLIMI.  THEH  TAKE  INVERSE  FTT  OF 
ROW  33  ONLY,  SINCE  IBIS  IS  THE  ROW  CONTAINING  CENTRAL 
INFCGMATIOH.  ROW  33  BECOMES  THE  NEXT  ROW  OF  THE  FINAL 
OOTFOT  MATRIX. 

DO  *13  COL  -  l.NCOLS 
DO  *05  ROW  -  l.NROWS 

CWK(ROW)  -  CONJG( WORK (COL .ROW) ) 

CONTINUE 

CALL  FFT2C  (CWK.M.IWK1) 

DO  A 10  ROW  -  l.NROWS 

WGRK(OOL.ROW)  -  CWK(ROW) 

CONTINUE 

CONTINUE 

DO  4*0  COL  -  l.NCOLS 

CWK(COL)  -  MGRK(COL.H) 

CONTINUE 

CALL  FFT2C  (CWK.M.IWK1) 

R123  -  NROWS  *  NCOLS 
C123  -  CMPLXCR123 ,0.0) 

DO  A 30  COL  -  l.NCOLS 

RESULT (COL , I )  -  ABS(REAL< <CONJG<CWK<COD )/C123 > ) ) 
CONTINUE 

CONTINUE 

ELSEIFtFTYPE  .EQ.  2)  THEN 
TEMPI  «  A  *  C  /  2 
TEMP2  -  A2  *  C2  /  A 
CALL  OOTOXm.l) 

WRITE<*,*)  ’CALCULATING  ROW  ’ 

DO  300  I  -  STEP+1, NROWS 
CALL  GOTOXYd,  18) 

WRITE(*. 101)  I 

DO  310  ROW  -  l.H 
DO  310  COL  -  l.H 

TEMP3  -  ZFRIME ( I )  *  SQRT(ABS(RH02(COL,ROW)**2-TEMP2) ) 
IF  ( ABS ( RH02 ( COL , RCW ) )  ,GT.  TIMP1)  THEN 
CALL  (MBS  JN  ( TEMPS  ,  N ,  TEMFA  ,  IER ) 

ELSE 

CALL  (MBSIN(TEMF3,N, TEMPA , IER) 

ENDIF 

WORK (COL, ROW)  -  TEMPA  *  EXP(-A*C2*TIME(I) )  * 

&  PULSE (COL, ROW) 

CONTINUE 


I 


J  -  0 


DO  320  SOW  -  H+l, SHOWS 
J  -  J  +  2 
DO  320  COL  -  l.H 

WORK (COL , HOW)  -  WORK(COL.ROW-J) 


320  CONTISUE 

C 

Q  A**************************************************************** 

C  FLIP  LEFT  SIDE  TO  RIGHT  SIDE 

C  A**************************************************************** 


J  -  0 

DO  313  COL  -  H+l.NCOLS 
J  -  J  +  2 

DO  313  ROW  -  1, SHOWS 

WORK (COL, ROW)  -  WORK ( COL -J. ROW) 

315  CONTINUE 

C 

C  *•*•******••********•**••*«,•****•******•********«.*•*•*••••*•*, 
C  TAKE  INVERSE  FFT  OF  EACH  COLIMf.  THEN  TAKE  INVERSE  FFT  OF 

C  ROW  33  ONLY,  SINCE  THIS  IS  THE  ROW  CONTAINING  CENTRAL 

C  INFORMATION.  ROW  33  BECOMES  THE  NEXT  ROW  OF  THE  FINAL 

C  OUTPUT  MATRIX. 

C 

DO  613  COL  -  l.NCOLS 
DO  603  ROW  -  1.NR0WS 

CWK(ROW)  -  CON JG( WORK (COL, ROW) ) 

605  CONTINUE 

CALL  FPT2C  (CWK,M,IWK1> 

DO  610  ROW  -  l.NROWS 

WORK ( COL , ROW)  -  CWK(ROW) 

610  CONTINUE 

613  CONTINUE 

C 

DO  6*0  COL  -  l.NCOLS 

CWK(COL)  -  WORK ( COL, H) 

6*0  CONTINUE 

CALL  FFT2C  (CWK.M.IWK1) 

C 

R123  -  NROWS  *  NCOLS 
C123  -  CMPLX(R123 ,0.0) 

DO  650  COL  -  l.NCOLS 

RESULT (COL, I)  -  ABS(REAL( (CONJG(CWKtCOL) )/C123) ) ) 

850  CONTINUE 

C 

300  CONTINUE 

C 


END  IP 


srr  bows  1  to  step  to  zero  to  simulate  step  function 
******************«*•**»— ****** «««««««  ****************  ««««««««  ** 
DO  500  SOW  -  1 , STEP 
DO  300  COL  -  l.NCOLS 
RESULT (COL, ROW)  -  0.0 

corn  hue 


OUTPUT  FINAL  ARRAY  TO  DISK 

a**************************************************************** 

DO  310  I  -  1.NR0WS 

WRITE (3, 102)  (RESULT (J, I) , J  -  l.NCOLS) 

CONTINUE 

FORMAT  (12) 

FORMAT  (64F16.7) 


CLOSE  REMAINING  OPEN  FILES 


CLOSE(2) 

CLOSE(3) 


REMIND  USER  OF  DATA  FILES  CREATED 


CALL  CLRSCR 
WRITE(*,*)  'FINISHED.' 

WRITE!*.*) 

WRITE!* , * )  '  THE  FOLLOWING  ASCII  TEXT  FILES  ARE  IN  THE  DEFAULT' 
WRITE!*,*)  '  DIRECTORY:' 

WRITE!* , *) 

WRITE!*, *>  '  FILENAME  DESCRIPTION’ 

WRITE(*,*)  '  . . . . . 

WRITE!*,*)  ’  ’, FRAME 1, '  INPUT  IMPULSE  FUNCTION’ 

WRITE(*, *)  ’  ' , FNAME2 , ’  SUhWARY  INFORMATION’ 

WRITE(* , *)  ’  ’ ,FNAME3 , ’  OUTPUT  ARRAY' 

WRITE(*,*) 

WRITE!*,*) 


*  SUBROUTINE  OOTOXY 


C 

C 

C 

C 

C 

C 


*  MOVES  CURSOR  TO  LINE  X,  COLUMN  Y  * 

*  CALL:  CALL  GOTOXY(X,Y)  » 

*  * 
*v<MHh**** **************** ***************************************** 


SUBROUTINE  GOTOKY(X,Y) 

CHARACTER* 1  C1,C2,CS,C8,LC(5) 

CHARACTER* 5  CBUFF 
INTE6ER*2  IC(4),L,X,Y 

EQUIVALENCE  (Cl, IC(1) ) . <C2, IC(2) ) . (C5, IC(3) ) . CC8 , IC(4) ) , 
+(CBUFF,LC(1)) 

DATA  IC/16#1B, 18#SB, 1S#3B, 16466/ 

C 

L  -  10000+100*X+Y 
C 

C  ***  WRITE  ESCAPE  CODES  TO  CHARACTER  BUFFER 
WRITE (CBUFF , 2)  L 
2  FORMAT(IS) 

C 

C  ***  WRITE  ESCAPE  CODES  TO  THE  DISPLAY 

WRITE(*,1)  C1,C2,LC(2),LC(3),C5,LC(4),LC(5),C8 
1  FORMAT ( IX, 8A1,\) 

C 

C  **«  END  OF  SUBROUTINE 
RETURN 
END 


C 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  CLRSCR  * 

SUBROUTINE  TO  CLEAR  THE  DISPLAY  * 

CALL:  CALL  CLRSCR  * 


SUBROUTINE  CLRSCR 
CHARACTER*!  C1,C2,C3,C4 
INTEGER*2  IC(4) 

EQUIVALENCE  (Cl, IC(1) ) . (C2, IC(2) ) , (C3 , IC(3) ) , (C4, IC(4 ) ) 
DATA  IC/18#1B,1643B.16#32,16#4A/ 

C 

C  ***  WRITE  ESCAPE  CODE  TO  DISPLAY 
WRITE(*,1)  C1,C2,C3,C4 
1  FORMAT ( IX, 4A1) 

C 

C  «**  END  OF  SUBROUTINE 


*  SUBROUTINE  TO  shift  a  squash  matrix  -  EXCHANGES  quadrants  ore  * 

*  AHD  THREE,  AHD  TWO  AMD  FOUR.  * 

*  USE:  CALL  SHIFT (X , BROWS , NCOLS )  * 


WHERE  X  IS  A  SQUARE  ARRAY  -  REAL  OR  COMPLEX  * 

HR0H5.RC0LS  IS  THE  ARRAY  SIZE  * 


SUBROUTINE  SHIFT  (X, NROWS, NCOLS) 


INTEGER  NROWS, NCOLS, ROW, COL, R2,C2 
COMPLEX  X(NCOLS, NROWS). T1.T2 

R2  -  NROWS/ 2 
C2  -  HCOLS/2 

DO  10  ROW  -  1,R2 
DO  10  COL  -  1.C2 
II  -  X(COL.ROW) 

X(COL.ROW)  -  X ( COL+C2 , ROW+R2 ) 
XCCOL+C2 , R0W+R2 )  -  T1  . 

T2  -  X ( COL+C2 , ROW) 
X(COL+C2,ROW)  -  X ( COL , ROW+R2 ) 
XtCOL .ROW+R2 )  -  T2 
CONTINUE 
RETURN 
END 
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